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One of the challenges with quantum simulation in ion traps is that the effective spin-spin exchange 
couplings are not uniform across the lattice. This can be particularly important in Penning trap 
realizations where the presence of an ellipsoidal boundary at the edge of the trap leads to dislocations 
in the crystal. By adding an additional anharmonic potential to better control interion spacing, and a 
triangular shaped rotating wall potential to reduce the appearance of dislocations, one can achieve 
better uniformity of the ionic positions. In this work, we calculate the axial phonon frequencies 
and the spin-spin interactions driven by a spin-dependent optical dipole force, and discuss what 
effects the more uniform ion spacing has on the spin simulation properties of Penning trap quantum 
simulators. Indeed, we find the spin-spin interactions behave more like a power law for a wide range 
of parameters. 


I. INTRODUCTION 

The idea of a quantum simulator, where a complex 
many-body quantum system is emulated in a controlled 
analog quantum computer and the results of the simula¬ 
tion are read off of the computer by measuring different 
properties as a function of time, originates with work 
from Richard Feynman [lj in the early 1980s. Cirac 
and Zoller [2\ showed how ion traps driven by a spin- 
dependent optical dipole force could realize quantum 
computers. Porras and Cirac Q further described how 
one could perform quantum simulations in a Penning 
trap. One of the issues with these quantum simulations is 
that the ions are not spaced uniformly. On the one hand, 
this leads to nonuniform effective spin-spin couplings be¬ 
tween the ions, on the other hand, in one-dimensional 
linear Paul traps, it leads to the linear to zig-zag transi¬ 
tion, which limits the number of ions that can be held in 
the trap in a one-dimensional linear configuration. It was 
quickly realized that by adding an anharmonic potential, 
which pushes together the farther out ions preferentially 
when compared to the central ions, one can achieve a 
more uniform arrangement, and the precise potential for 
perfectly uniform trapping is known for the linear Paul 
trap ( 5 } . Surprisingly, one can achieve quite uniform crys¬ 
tals by just adding a quartic potential on top of the con¬ 
ventional quadratic trapping potential. This ideology has 
been extended to the Penning trap by Dubin i5j, where 
he also included a triangular-shaped rotating wall poten¬ 
tial to reduce dislocation formation, which occurs in two- 
dimensional Penning traps when the boundary potential 
does not have the same symmetry of the underlying ionic 
lattice. 

In this paper, we extend the analysis of Dubin to deter¬ 
mine the behavior of different numbers of trapped ions, 
different wall potentials and different rotation rates, to 
determine the stability of these ionic crystals. We fur¬ 


ther calculate the axial phonon modes and from them 
the effective spin-spin interactions induced by a state- 
dependent optical dipole force. We end by discussing the 
feasibility of the triangular wall for quantum simulation 
with ionic crystals in the Penning trap. The organization 
of this paper is as follows: in Sec. II, we describe the the¬ 
oretical background for the calculations. In Sec. Ill, we 
present the numerical results for the calculations of the 
spin-spin couplings of the ions. In Sec. IV, we present 
our conclusions. 


II. THEORETICAL FORMULATION 

The Penning trap confines ions by using an electro¬ 
static potential that pushes the ions towards the plane 
with z = 0, and also pushes the ions outwards, radially. 
A large static magnetic field curves the radial motion 
into circles, which result in a trapped ion crystal (after 
taking into account the Coulomb repulsion of each ion). 
An additional rotating wall potential, with specified an¬ 
gular symmetry, is then applied to control the shape of 
the crystal and its rotation rate. While many different 
ions can be employed in a Penning trap, we will focus 
on the realization with two liyperfine levels of 9 Be + ions, 
\ 2 S 1 / 2 ,mj = 1/2) and \ 2 S 1 / 2 ,mj = —1/2), localized to a 
single plane. An extensive description of such a set-up 
can be found elsewhere (f|. Cold atoms condensing in 
such crystals are candidates for building quantum simu¬ 
lators, owing to the ease with which these systems can be 
prepared for large number of ions, and the precise quan¬ 
tum control of individual ions these systems afford !7]. 

In actual experiments, the ionic crystal often acquires 
an increasing number of impurities as a function of time, 
such as BeH + that form due to collisions of the beryllium 
ions with hydrogen molecules that are in the dilute back¬ 
ground gas. While a theoretical treatment that includes 
the effect of these impurities is possible [§], we consider 
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only the clean limit here, where there are no impurities. 
This simplifies the analysis below. In addition, exper¬ 
imental protocols to purify the systems may make this 
effect less important Q. 

The theoretical treatment of the equilibrium positions 
and normal modes of the Penning trap requires a careful 
analysis employing standard classical mechanics and then 
an appropriate quantization scheme^ Details for how to 
do this have already appeared EMI- Here, we provide 
a quick summary of that formalism to establish our nota¬ 
tion and to show how the approach needs to be modified 
for the more uniform triangular crystals that can be gen¬ 
erated in an extra anharmonic potential with a rotating 
triangular wall. We begin with the ion Lagrangian in the 
laboratory reference frame which satisfies 


N 


j =1 



+ eA-Tj 


(1) 


where N is the total number of ions, e is the (positive) 
unit charge of an electron and m is the mass of a 9 Be + 
ion. The symbol r j = ( Xj,yj,Zj) is the position vector 
for the jth ion in Cartesian coordinates, is the total 
scalar potential acting on the jth ion (including the ro¬ 
tating wall potential), and Aj = (B x ij)/2 is the vector 
potential in the symmetric gauge for the uniform axial 
magnetic field B = B z z with ( B z > 0). The scalar po¬ 
tential 4>j(t) includes the potentials that trap the ions 
and the mutual Coulomb repulsion between the ions. It 
can be expressed as follows: 
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where ui e f f = y a; c S2 — fl 2 — is the effective trap¬ 
ping frequency in the rotating frame of the crystal with 
w c = eB z /m (see below), Vo is the amplitude of the static 
quadratic potential from the Penning trap electrodes, C 4 
is the strength of the additional fourth-order anharmonic 
trapping potential (also coming from the Penning trap 
electrodes), V\y is the amplitude of the triangular rotat¬ 
ing wall potential, D > 0 is the rotating wall angular 
frequency (which rotates about the 2 -axis), and k e is the 
Coulomb force constant. Here, r k j = |rk — rj | is the 
interparticle distance between the fcth and jth ion and 
is given by y/(x k - Xj) 2 + (y k - %j 2 + (z k - Zj) 2 , pj is 
the polar coordinate radius for the jth ion and is given 

by Pj = \J'%] + Uj , and 0j is the polar angular coordi¬ 
nate for the jth ion which is given by 0j = tan -1 ^- /xj). 
The rotating wall potential makes the potential (pj time- 
dependent in the laboratory frame. The solution for the 
ion positions is simplified by transforming to the equiv¬ 
alent equilibrium problem in the rotating frame with 
angular speed (where the effective trapping potential 


becomes time-independent). Transforming to the rotat¬ 
ing frame, we arrive at the following time-independent 
rotating-frame potential for the jth ion (which is con¬ 
fined to the plane with z = 0): 

e( t>j = \mul f f\(pf) 2 + C 4 (pf) 4 ] 

+ Vw[(xf) 3 -3xf(yf) 2 ] + ^-J2^R ( 3 ) 

fe/j T jk 

where rf = (xf, yf, zf) is the transformed set of coor¬ 
dinates for the rotating frame. The fl-dependence of the 
effective trapping frequency is due to the way velocities 
transform for rotating frames, which now includes the ef¬ 
fects of potential terms from the centrifugal and Lorentz 
forces as well. 

Comparing Eq. 0 with a similar expression appearing 
in Ref. [12, we see the following differences: (i) there is 
an additional anharmonic, fourth-order term that causes 
the ions in the outer regions of the crystal to be pushed 
in more strongly, and hence counteracts some of the in¬ 
homogeneities that occur due to increasing interion dis¬ 
tances as we move outwards; (ii) the angular shape of 
the rotating wall term has now been adjusted to an 
l = 3 angular harmonic as the crystal condenses into 
a triangular lattice and there is less frustration at the 
edges if the rotating wall has a symmetry that matches 
that of the underlying crystal facets Q. We note that 
the additional/modified terms retain their form under a 
transformation from the laboratory frame to the rotating 
frame. We also remark that the unconventional choice of 
(with the coefficient dependent on f l) for 
the anharmonic term was made in anticipation of the 
simpler and rather standard form we get in Eq. ©. 

The ions in a Penning trap crystal do not always crys¬ 
tallize in a two-dimensional plane. The following approx¬ 
imate criterion is usually required to be satisfied for a 
planar only configuration: 


2eV 0 

(eB z Q — mfl 2 — eVo ) 
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( 4 ) 


which basically ensures that the restoring force in the 
axial direction is several orders of magnitude bigger than 
in the radial direction, and the crystal lies in the z = 0 
plane only. 

The triangular rotating wall also introduces an 
anisotropy in the radial potential, with deconfinement 
along certain directions which, if strong enough, can lead 
to particle loss. For the more familiar case of l = 2, 
this is reflected in the increasing eccentricity of the ellip¬ 
tical equilibrium structures that suggests deconfinement 
along the ’’weak” axis. This is expressed in terms of an 
approximate criteria in Ref. I12j for the simplest rotat¬ 
ing wall. For l = 3, this becomes difficult to express 
in terms of a single criterion and we resort to a numer¬ 
ical calculation to find the triangular ” separatrix”, the 
contour lines of the radial part of the gradient of the 
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potential function in Eq. ©, with the Coulomb repul¬ 
sion terms ignored. These are shown for varying ratios 
of strengths of the rotating wall term to the effective ra¬ 
dial confinement strength, Vw/w^ff in Fig. [T] We notice 
that unlike the quadrupole rotating wall potential, decon¬ 
finement happens for the triangular wall at large enough 
distances for all wall amplitude strengths. We see that 
as the strength of the triangular rotating wall increases, 
the separatrix moves closer to the center of symmetry of 
the crystal, with apparent deconfinement centered along 
the 0 = 0, 27r/3 and 47r/3 axes. This has a noticeable 
effect on the ’’shape” of the Penning-trap crystals, which 
reduces the dislocations in the crystal and helps maintain 
the uniformity. Of course, a full description of stability of 
these planar structures requires inclusion of the Coulomb 
terms. The accurate quantitative description of stability 
requires solving for the equilibrium positions presuppos¬ 
ing a planar arrangement of ions, and then showing that 
the eigenvalues of the normal modes of oscillation about 
these equilibrium positions are all positive; that is, the 
phonon normal mode frequencies are all real. 

The full, transformed Lagrangian for the rotating 
frame is 


N r 

j =1 


2 m|r J 


R| 
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(5) 

where B e ^ (Q) = B z — 2VLm/e is the fl-dependent effec¬ 
tive magnetic field in the rotating frame. The modifica¬ 
tion of the magnetic field is due to velocity dependent 
terms in the laboratory frame Lagrangian. This affects 
the oscillating normal modes of the planar motion when 
the ions are far from their equilibrium positions. How¬ 
ever, as we will see soon, this does not have an effect on 
the axial modes. This observation greatly simplifies the 
normal-mode analysis for the axial modes. 

To find the stable spatial configuration of the ions, we 
minimize the effective potential energy in the rotating 
frame of reference. This is a challenging optimization 
problem to solve for in two (and higher dimensions), es¬ 
pecially since different configurations, separated by large 
potential barriers, can have local minima in the poten¬ 
tial energy function with small energy differences to the 
global minimum. We follow the previous treatments of 
this problem, where the experimental indication of the 
fact that the ions condense in a triangular lattice in a 
single plane is used to construct the optimized solution 
that lies close to a perfect triangular lattice. 

We construct an initial, trial solution based on the 
“closed-shell” approximation as in Ref. [TH but with the 
important difference that the overall “shape” is triangu¬ 
lar and not hexagonal, as it was in the previous solu¬ 
tions. This is to reflect the fact that the overall shape 
of the crystal is dictated by the equipotential lines of the 
rotating wall term, which in this case is triangular. 

We then proceed to calculate the collective normal 
mode excitations of the crystal. The ion Lagrangian (in 
the rotating frame, all R superscripts are dropped for 


clarity) is expanded via a Taylor series about the previ¬ 
ously calculated equilibrium positions of the ions up to 
quadratic order. The ion coordinates are, for the purpose 
of the expansion, written as r, (t) = R? + 5Rj(t), while 

for the ion velocities we write rj ( t ) = <5Rj ( t ). Because we 
are expanding about equilibrium, we can drop the linear 
terms in the coordinates, and we find 
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where the Cq is due to the equilibrium state and the 
quadratic terms are due to fluctuations away from equi¬ 
librium, which we henceforth call C p h for the phonon 
Lagrangian. The Lorentz force due to the external mag¬ 
netic field lies in the xy plane, and the potential energy cf>j 
is clearly seen to be separable in cylindrical coordinates. 
This means that the axial phonon Lagrangian can be de¬ 
coupled from the planar phonon Lagrangian, and there 
is no harmonic coupling between the planar and axial de¬ 
grees of freedom. Therefore, we can study the axial and 
planar modes independently ( C p h = + C pl £ nar and 

we can solve just the equations of motion for the axial or 
the planar modes independent of the other). 

We examine only the axial modes in this work. This 
is due to the fact that the planar modes have a complex 
structure owing to a coupling of the ion motion in the 
x and y directions, the appearance of velocity-dependent 
forces, as well as complexities introduced by rotation of 
the ion crystals as observed in the laboratory frame. The 
fact that there is no harmonic coupling between axial and 
planar directions of the crystal allows us to exclude the 
planar modes from our discussion henceforth. Restricting 
to the axial modes only, is further supported by the fact 
that the simplest form of quantum simulation works on 
driving the axial modes with a state-dependent optical 
dipole force. 

The axial Lagrangian is then given by 


r axial 
^ph 


1 

2 


E 771 ( SJ *k) 
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k =1 


\ E K ^SR z k 
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(7) 


where the spring constants satisfy 


d 2 C 

K zz = - 
3k dR~dRl 


( 8 ) 


The absence of cross terms in the velocity part of the 
Lagrangian can be easily seen from the sum-of-squares 
structure of the kinetic energy along the 2 -direction and 
the fact that there are no velocity-dependent forces in 
the ^-direction. 

An explicit calculation for the matrix elements of the 
K zz gives the following: 


K * = 


2eVo k e e*'£ k ', k >& (E,) 3 
7 „ r 2 1 
e WJkF 
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j ^ k 
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where R°- k = |R° — R°| is the distance between ions 
located at their respective equilibrium positions in the 
rotating frame. We see that the axial stiffness matrix 
is Hermitian and symmetric, and is independent of the 
anliarnronic or wall potentials. 

To solve for the axial ion normal modes, we apply 
the Euler-Lagrange equations to the axial phonon La- 
grangian in Eq. 0: 

N 

rnSRj + y K jk^ Z k = 0 , j = l,2,....N. ( 10 ) 

k =1 

which, on substitution of the eigenvector solution ansatz 
SRj(t) = bj 1 ' cos[(jj Z i,(t - t 0 )], gives 

N 

E [mu 2 zv S jk - K$\ tif? = 0 , j,u = 1,2 ,...N ( 11 ) 

fc=l 

where uj zv is the normal mode eigenfrequency and b z k is 
the i/th axial normal-mode eigenvector. The eigenvalue 
problem is quadratic, but we can easily map it onto a 
linear eigenvalue problem by setting the eigenvalue ac¬ 
cording to \ zv = mu>1 v . We can then solve for the eigen¬ 
values and eigenvectors numerically in MATLAB. The 
eigenvectors b zu are real, N-tuples whose norm has been 
set to unity by convention. The eigenvalues \ zv are pos¬ 
itive for stable normal modes and negative for unstable 
normal modes. 

The quantization of the normal modes is completely 
standard: we first identify the positions Q v and mo¬ 
menta P v associated with each phonon mode as canon¬ 
ically conjugate, and promote the relation given by the 
Poisson bracket {Qv, Pv'} = Sw' to the commutation re¬ 
lation for the operators Qv and P v >, [Q v ,Pv'] = ihSw’- 
To calculate the canonically conjugate variables for the 
phonon modes, we make the transformation SR z (t) = 
Y^v where £„ are the normal coordinates for each 

phonon mode v. We see that the Lagrangian assumes the 
following diagonal form: 

1 N 

£ a p t al = £ <*&£)• ( 12 ) 


Hence, we calculate the conjugate momenta as follows: 


ftr axial 
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The Hamiltonian is then expressed as 
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To quantize the normal modes, we identify that the 
Hamiltonian H ,z k lal is a sum of simple harmonic modes 
with frequencies ui zv . We now introduce creation and 
annihilation operators as follows: 


muj zv 
2 h 
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-P 


axial 


and 


mu zv 
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muj zv 
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mcj z 




(15) 
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Hence, the quantized Hamiltonian operator can be ex¬ 
pressed as 
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J j 

tL, 


axial 

ph 


y hu) z 


i 


zv l ,b zv 


(17) 


where n zv = a' zv a Z v is the number operator. The op¬ 
erator for displacement along the ^-direction can be ex¬ 
pressed in terms of the creation and annihilation opera¬ 
tors as follows: 


6R J = j2 b r 

v=\ 


2muj z 


[a\v + a z 


(18) 


We also note that the form of the Hamiltonian derived for 
the axial modes here is invariant when we transform co¬ 
ordinates from the lab to the rotating frame, since the ion 
oscillations are only along the rotation axis (^-direction) 
and hence are not influenced by rotation of the coordi¬ 
nate axes. 

We now need to calculate the effective spin-spin cou¬ 
pling between the ions, generated by the spin-dependent 
optical dipole force. This analysis has been done in detail 
elsewhere [1, HH, HH , which we utilize here. The effective 
spin Hamiltonian is dictated by a time-dependent Ising 
spin Hamiltonian, 
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u(t) = y . 

j,r =i 


(fi' (t) a j a j ’> 


(19) 


where the Ising spin-spin coupling between sites j and j' 
is given by 
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( 20 ) 


Here Fo is the magnitude of the optical dipole force, 
and fi is the beat-note frequency corresponding to the 
frequency difference of the two off-resonant laser beams 
being applied to the trapped ion crystal. We see that this 
expression relates the strength of the Ising-like coupling 
between ions to the phonon mode properties ui zv , b zu and 
b z Y , which are calculated from the classical, normal-mode 
analysis described above. The time-averaged spin-spin 
couplings are given by the first term in Eq. We 

can think of the effective spin-spin Hamiltonian as the 
expression in the parenthesis of Eq. m with the time- 
dependent spin-spin interactions replaced by the time- 
averaged ones. 
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III. NUMERICAL RESULTS AND ANALYSIS 


We consider 9 Be + ions localized in a plane by the Pen¬ 
ning trap potential defined in Eq. ©. We character¬ 
ize the strength of the end-cap potentials Vq that af¬ 
fect the axial trapping by a characteristic angular fre¬ 
quency w z , such that eVo = ^mu> z . This is fixed for 
the purpose of all our numerical calculations at the value 
w z = 2tt x 795 kHz, a typical value used in experiments. 
We normalize subsequent frequencies in terms of lo z . The 
experiments at NIST typically run at rotational frequen¬ 
cies fl = 0.0579w z , and we have concentrated on regions 
close to this value in our calculations for experimental 
relevance. The cyclotron frequency w c associated with 
the magnetic field is defined as w c = eB z /m. Fixing 
B z = 4.5 T, we get w c = 9.645w z . The beryllium atom 
has an atomic mass m = 9.012182 a.u. and a positive 
unit charge e = 1.60217646 x 10 -19 C. 

For the strength of the anliarmonic term, we use the 
value in Ref. [1, where C 4 = 1 in the following form of the 
potential s: 
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2 muj e ff 


pi 


3 ^ Pi 

-L-4-y 
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Pi cos(3 0i) 
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( 21 ) 


where r p is the plasma radius parameter whose value is 
taken to be r p = 0.01049 cm. We now define a typical 
length and energy scale, 


lo = 



and E = 


( 22 ) 


Henceforth, we express all lengths in units of Iq and all 
energies in units of E. The potential then becomes 
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= - = y 

E ^ 

1=1 

+ Vw(Xi - 3 xml) 


Iff {pi + c^i) + -z X/ 


2 *-?. ri 


(23) 


where the dimensionless parameters C 4 and Vw are 
given by the relations C 4 = 31qC 4/(8r p ) and Vw = 
HvW( mtJ z r p)- All lengths (Iq), energies (E) and fre¬ 
quencies (u> z ) appearing in Eq. E3l) are expressed in di¬ 
mensionless units. For the choice of C 4 = 1, we have 
C *4 = 0.002472. For the strength of the wall potential, we 
use two different values, Vw = 0.0025 and Vyy = 0.0040. 
These have been chosen to show clearly the effect of vari¬ 
ation in the wall strength on key ion-crystal characteris¬ 
tics, while also ensuring that all crystal structures have 
a stable equilibrium. 

We report the rotational frequency Q in terms of 
oj e ff = — fl 2 — 1 / 2 , normalized with respect to 

oj z . We stick to regions close to the experimental value 


of 11 = 0.0579w z , which translates as u> e ff = 0.2339w z . 
We use values for w e // in the range 0.21 u z — 0.25w z in 
our analysis. 

It is also useful to consider the stability of the crys¬ 
tal under the trap potentials we have used. Because the 
rotating wall potential varies as the third power of the 
coordinates, there is no deconfinement frequency as in 
the case of the quadrupole wall, l = 2 trap. To probe the 
stability under deconfining forces, we look at the radial 
component of the force on the «th ion due to the trap po¬ 
tential, excluding the Coulomb potential, which is given 

by 


F r - -uil ff 


Pi 


2 C 4 3 
~2 Pi 

U eff 


3Vw 
W Iff Pi 


(x 3 i - 3 Xiyi) 


(24) 

If we plot the locus of points where this function 
becomes zero for various values of Vw, we see that 
we get three regions, each centered along the 8 i = 
— 7 r/ 3 , 7 r /3 and 7 r axes. For increasing strength of the 
rotating wall potential, they move closer to the origin 
and the radius of the separatrix (the smallest distance to 
these unstable zones along any axis) is seen to decrease. 
Hence, the deconfinement increases for very high values 
of the rotating wall potential, which is what we expect. 
In our analysis, we stay in regions where the extent of 
crystal is much smaller than the separatrix radius. 



FIG. 1. (Color online.) Contour plot of the magnitude of 
the radial restoring force F r = 0, for parameterized val¬ 
ues of the strength of the rotating wall potential (given by 
3 Vw/^eff = 0.508, 0.608, or 0.908). The values of this 
strength have been indicated by appropriately colored labels 
near each curve. The separatrix radius is roughly the short¬ 
est distance to the contour, from the center of the graph and 
we see that the separatrix radius decreases as we increase the 
strength of the rotating wall term. 

We also need to exclude unstable equilibrium configu- 
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rations (which are indicated by the nonpositivity of the 
eigenvalues \ ZIJ of the stiffness matrices of the axial vi¬ 
brations). We examine this in detail below, when we 
discuss our results on the axial phonon modes. 

To find the equilibrium configurations, we need to 
minimize the Hamiltonian of the crystal in the rotating 
frame, which boils down to finding the best minimum 
of the transformed potential function of Eq. (l23l) near a 
triangular lattice. We only concentrate on the solution 
we obtain starting from the closed-shell construction. In 
this case, we start with a seed lattice where we arrange 
ions in closed, triangular shells while also respecting the 
triangular lattice symmetry. The shell is closed if we can 
put all the ions in these complete shells. We relate the 
number of shells S to the number of total ions in the 
crystal N as 


S = 


2(N — 1) 



1 

2 


(25) 


If we cannot put the N ions in an integer number of 
shells, we put the outermost ions in an incomplete tri¬ 
angular ring according to the minimal potential energy 
at each of the outer ring sites. For the purpose of our 
discussion here, we pick N = 85 and hence S = 7. We 
arrive at the minima guaranteed under such a considera¬ 
tion of the seed lattice using a trust-region algorithm of 
the MATLAB Optimization Toolbox. The minimization 
procedure requires us to specify a locally calculated gra¬ 
dient of the potential, which can be input analytically by 
taking derivatives of the potential. The procedure iter¬ 
ates the minimization steps until the local minimum is 
found. 

The equilibrium configurations we obtain from such a 
procedure behave as we might expect (see Fig. [2]). We 
obtain structures that form a nearly perfect triangular 
lattice close to the center, and smoothly transition to the 
shape of the contour lines of the effective potential as we 
move radially outward to the edges. The edge effects 
cause the interionic distances to change as we move out¬ 
wards. Our strategy to counter these effects are two-fold, 
as discussed earlier: (i) We introduce a weak anliarmonic 
term whose strength is characterized by C 4 and (ii) we 
match the symmetry of the rotating wall with that of the 
condensed crystal. The first effect has been incorporated 
and fixed at a particular value, as discussed earlier. How¬ 
ever, we keep the strength of the rotating wall variable 
as it is only for a certain range of values of the strength 
of the rotating wall potential that the contour lines of 
the effective potential are triangular, and hence become 
least likely to cause edge distortions of the ionic crystal. 

This is clearly seen in the Fig. [2] where we show the 
progression of these structures for different values of the 
effective radial trapping strength u> e /f, for = 0.0025 
and Vyy = 0.0040. For the low-strength rotating wall, 
the structures (a)-(d) are nearly triangular and uniform, 
whereas the higher value of the rotating wall strength cor¬ 
responds to the more distorted structures of (e)-(h). For 


a fixed number of ions and fixed C 4 , we see the desta¬ 
bilizing effect of the decrease in separatrix radius with 
increasing rotating wall strength clearly on these crystal 
structures. We note here that in the limit of vanishing 
wall strength or high effective radial trapping frequency, 
the structures become more isotropic, with uniformly de¬ 
creasing nearest-neighbor distances as we move towards 
the edges. The other extreme (very small radial trap¬ 
ping frequency or high strength of rotating wall), the 
equilibrium configuration we obtain from a closed-shell 
construction shows that the ions are reduced to (three) 
pockets of stability and the structure is no longer closed. 
A detailed normal-mode analysis (see below) shows that 
these structures are in unstable equilibrium, and hence 
we can discard them. Another important observation is 
that a given progression of structures (for differing values 
of the radial trapping strength and increasing values of 
Vw) displays similar structures to those found at smaller 
trapping strengths, for higher values of the rotating wall 
strength. This fact will be important to arrive at struc¬ 
tures that show the maximum uniformity and also exhibit 
a stable equilibrium. 

Note that in these calculations, we fix the wall po¬ 
tential and then vary w e //. In doing so, we find that 
we are limited by how many ions we can hold in stable 
equilibrium. Because the radial deconfinement decreases 
as the rotational frequency w e // decreases, the effect of 
the rotating wall will become stronger if it also remains 
fixed. We have done this here to reduce the parameter 
space we explore. But, if one wants to examine larger 
crystals, then one needs to carefully tune the rotating 
wall strength as the rotational frequency is changed, as 
well as the strength of the quartic potential, to be able 
to continue to maintain stable equilibrium. These issues 
are discussed in Ref. @ 



FIG. 2. Equilibrium structures found for varying w e // and 
Vw- Panels (a)-(d) show the succession of crystal structures 
obtained by increasing w e // with Vw = 0.0025w*. Panels (e)- 
(h) present the structures for higher Vw = 0.0040oj^ for the 
same corresponding values of ui e f /. 


The equilibrium structures we obtain are markedly 
uniform, and this is born out in Fig. [3] where we plot 
the distance to the nearest neighbor for each ion in the 
crystal as a function of the central ions’ distance from 
the trap symmetry axis. These results are for the most 
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uniform stable structures we could obtain for both the 
l = 2 and l = 3 rotating walls. The nearest-neighbor 
distances have been calculated based on the Delaunay 
triangulation algorithm. Each point in the figure repre¬ 
sents the distance of the ion in question to an ion in the 
first nearest-neighbor shell. We focus on the first circle 
of nearest neighbors only. The values of the relevant pa¬ 
rameters, for both the l = 2 (quadrupole) and the l = 3 
(triangular) wall crystals, are indicated in the caption to 
Fig. [3] The larger spread of values for the quadrupole 
wall tells us that, on average, the triangular wall crys¬ 
tal is indeed more uniform spatially than the quadrupole 
wall crystal. Note that there are the same number of 
blue squares and red circles in Fig. [3] the uniformity of 
the triangular lattice has many of the symbols overlap. 
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FIG. 3. (Color online.) Nearest-neighbor distances d versus 
the distance of the ion from the trap symmetry axis p , where 
all distances have been normalized against the length scale 
lo- The red circles plot this for the Z = 3 triangular wall with 
an additional anharmonic potential for Vw = 0.0025oif and 
w e // = 0.25w z [equilibrium positions shown in panel (b)], 
while the blue squares represent this for the original l = 2 
rotating wall with just a harmonic potential for w e // = 0.06w z 
[equilibrium positions shown in panel (c)[. The parameters for 
the former trap are chosen so as to minimize the variance in 
nearest-neighbor distances, while a moderate wall potential 
is chosen for the latter for a valid comparison. The variance 
in the triangular wall lattice (b) is much smaller than in the 
quadrupole wall lattice (c). 


We next discuss the features of the normal modes of 
small oscillations of the ions. We first examine the pos¬ 
itivity of the eigenvalue spectrum of the stiffness matrix 
of axial vibrations K zz as a function of effective trap¬ 
ping strength, for high and low rotating wall potential 
strengths, in Fig. 2) We see that the eigenvalue spec¬ 
trum is real and positive for only a “band” of values 
of w e // (indicated in blue in the figure), and this band 
shrinks for the higher amplitude rotating wall potential. 
The positivity of all the eigenvalues indicates stability of 
the corresponding structures. 

This behavior is fundamentally different from that seen 
in the usual Penning trap crystals where we do not see a 
lower bound to the strength of the radial trapping, char¬ 
acterized by w e //, where the structures become unstable, 
although there is an upper limit. Hence, we can only talk 
of stable structures for certain ranges of w e // for the l = 3 
rotating wall (with an additional quartic potential), and 
this range gets narrower for increasing wall strength. A 
similar effect is seen if the number of ions in the trap 
is increased, in that the band of stability shrinks as we 
increase the number of ions in the trap. The value of 
N = 85 was the maximum number of ions we found that 
could be trapped with reasonably large bands of stabil¬ 
ity for the particular ranges of parameters that we chose. 
Note that more ions can be trapped by carefully choosing 
the rotation frequency w e //, the rotating wall potential 
amplitude, and the strength of the quartic potential, but 
we do not discuss these cases in detail here. 

Next, we plot the numerically obtained eigenfrequen- 
cies against their mode numbers, for crystal structures 
corresponding to two values of (stable) w e //, for both 
high and low rotating wall potentials. Roughly, we can 
see a trend similar to that of the quadrupole-wall crystal 
normal modes, where the primary dependence is on w e // 
and not on Vw- However, there are important distinc¬ 
tions to be made. The “band” structure of the eigenvalue 
spectrum causes the structure of these curves to change, 
and quite significantly at that, when we vary the rotat¬ 
ing wall from low to high strength. We see this in Fig. [U 
where for ui e ff = 0.20w~, the lower edge of the instability 
of the band shown in FigQ] shifts to the right when the 
strength of the rotating wall is increased, and this causes 
the eigenfrequencies for the higher wall strength to drop 
abruptly to values very close to zero. In this fashion, we 
see that the dependence on w e // is now superimposed on 
a dependence on the strength of the rotating wall due to 
a Vvy-variable bandwidth. 

The highest axial mode has a universal eigenfrequency 
equal to the angular frequency of the trapping strength, 
ui z for all values of w e // and Vw, and we see that all 
the branches converge to this point. This behavior is 
identical to that of the quadrupole-wall rotating crys¬ 
tal. The corresponding eigenmode is the well-known 
center-of-mass mode, where all ions move equal displace¬ 
ments that are in-phase with each other. This is because 
the center-of-mass motion does not cost any additional 
Coulomb energy and all ions have the same axial trap- 
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FIG. 4. (Color online.) Variation of the axial normal mode 
frequencies (in units of uj z ) with cf° r two wall strengths: 
(a) V w = 0.0025w z and (b) Vw= 0.0040^. 


ping energy, and hence their motion is independent of 
the strength of the rotating wall potential applied in the 
crystal plane. Also, other axial phonon modes will have 
frequencies lower than the center-of-mass mode, as the 
average distance between the ions increases when the 
wave vector is nonzero and there is a reduction in en¬ 
ergy due to the Coulomb repulsion. This is also why the 
axial eigenfrequency branches for higher co e f / he roughly 
lower in Fig. [5] 

Next, we discuss the eigenvectors of the axial mode 
phonons, which we can calculate immediately from the 
diagonalization of the stiffness matrix that yields the 
eigenfrequencies. Modes close to the center-of-mass mode 
are collective, where ions move with a long-wavelength. 
It is these long-wavelength modes that are important for 
purposes of quantum simulation, and we concentrate on 
these in our discussion here. In Fig.[6l we show maximum 
axial displacements of the ions, corresponding to the 
highest three axial modes, with the displacements nor¬ 
malized to unity and color-coded as indicated in the adja¬ 
cent colorbar. We do this for the two different strengths 
of the radial trapping strength w e //, for both the high 
and the low rotating wall potentials Vw- The highest ax- 



FIG. 5. (Color online.) Eigenfrequencies of the axial phonon 
modes. The green and red symbols represent data for a 
strong rotating wall, Vw = 0.0040u; z , and a weak wall 
Vw = 0.0025tJ^, respectively. The values of the corresponding 
effective trapping frequency io e ff are indicated by the labels 
near each curve (hollow symbols, u) e ff = 0.20oj z and solid 
symbols, u} e ff = 0.26u; z ). Note how the strong rotating wall 
and high trapping frequency case is nearly unstable. 


ial mode corresponds to the center-of-mass motion. The 
other two modes (so-called tilt modes) are nearly degen¬ 
erate, and are seen to have similar nature even when 
anisotropy due to the rotating wall is dominant. This is 
also seen in Fig. [5J where we see that the two modes just 
below the center-of-mass mode have the same eigenfre¬ 
quencies. This behavior of the penultimate axial modes 
is different from the quadrupole-wall crystal, where this 
degeneracy of the modes is lifted under the correspond¬ 
ing anisotropic rotating wall, and an additional mode is 
sometimes introduced between them. 

Finally, we examine the strength of the effective spin- 
spin coupling that results from applying a spin-dependent 
dipole optical force detuned close to the axial phonon 
modes. We focus only on detuning frequencies <5 = fi— uj z 
to the blue of the center-of-mass mode (6 > 0) where 
we expect to find a power law dependence of the spin 
exchange on the distance between the spins in the lat¬ 
tice jUJ. This behavior has been predicted and verified 
in experiments, and we expect it to hold even for this 
triangular-wall anharmonic crystal. In fact, the objec¬ 
tive for making the lattice distances more uniform, is 
to also make the spin-spin couplings between adjacent 
spins more uniform, and this is would follow immedi¬ 
ately if there exists a power law relation between the two 

(Jjj’- J 0 /|R2-R$,|“)._ 

The parameters of axial phonon modes discussed al¬ 
ready are used to calculate the static spin-spin interac¬ 
tion Jjj* between the spins of ions j and j', based on 
Eq. (EU1) . In Fig. [3 we plot this static interaction strength 
(expressed on a logarithmic scale) as a function of the 
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FIG. 6. (Color online.) Three highest-frequency axial eigen¬ 
vectors for various trapping strengths 0J e ff and rotating wall 
potentials Vw- Panels (a) — (/) correspond to a low ro¬ 
tating wall strength Vw = 0.0025^, (a) — (c) represent¬ 
ing eigenmodes of the crystal structure corresponding to 
w e // = 0.21w z and ( d ) — (/) representing the crystal structure 
for w e // = 0.24 uj z . Panels ( g ) — (Z) correspond to a stronger 
rotating wall strength Vw = 0.0040w^, ( g ) — ( i ) corresponding 
to w e // = 0.21w z and ( j)—(l ) corresponding to w e // = 0.24w z . 


distance between the ions, for various values of detun¬ 
ing /i larger than the center-of-mass frequency. We do 
this for the stable crystal structure exhibiting minimum 
variance in the nearest-neighbor distances. We see a be¬ 
havior very similar to the quadrupole-wall potential. A 
uniform J ^, independent of indicates that the detun¬ 
ing laser excites only the uniform center-of-mass mode. 
As we move away from the center-of-mass mode, an in¬ 
creasingly large number of eigenmodes participate in the 
coupling, and we see a clear trend in the values, such 
that Jij oc In the limit of large detuning, we have 
dipole-dipole interactions where a tends to a value of 3. 
For small detunings, we have the all-to-all case of a —» 0. 
There are small departures from the power law behavior 
for intermediate detunings, while the very small and very 
large values of detuning show excellent agreement with 
the power law. Note that these results are indeed much 
more uniform than what was found for a a quadrupole 
rotating wall @. 

In Fig. [5J we plot the fitted power-law exponent a ver¬ 
sus the strength of the detuning away from the center- 
of-mass mode. The trend we see here is similar to the 
one seen in calculations for the quadrupole-wall poten¬ 
tial. We see a faster approach to the dipole-dipole limit 



FIG. 7. (Color online.) Time-averaged spin-spin coupling co¬ 
efficient Jij/J ( J = FQ/muOz) versus the distance between 
the ions rij/lo = |R° — Rj|/Zo on a log-log plot. We plot 
the stable crystal structure exhibiting minimum variance in 
nearest-neighbor distances, corresponding to Vw = 0.0025tu z 
and ijJe/f = 0.25oj z ■ The power law exponent a and the 
strength of the detuning away from the center-of-mass mode 
J are both indicated near each curve. 



FIG. 8. (Color online.) Fitted exponent of the power law a (of 
the Jij coefficients as a functions of the distance rij ) plotted 
against the strength of detuning away from the center-of-mass 
frequency, 8, for the same set of trap parameters and notation 
as in Fig. [5] 


(a = 3) for both smaller effective (radial) trapping fre¬ 
quencies, and weaker strengths of the triangular rotating 
wall potential, just like for the quadrupole rotating wall. 

We have already noted that there are deviations from 
the power law behavior for intermediate values of detun¬ 
ing, and this is apparent in the spread of values away 
from the linear fit in Fig. [7] (especially for <5 = 10 -1 u; z ). 
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FIG. 9. Normalized oot mean square deviation for the fits 
of the spin coupling constants to a power law for different 
detunings to the blue of u z . 


To explore these deviations in more detail, we plot the 
normalized root mean square deviation (RMSD), defined 

by, 


the limit of dipole-dipole interactions, we get a value of 
a very close to 3, with the normalized RMSD close to 0. 

We note that such behavior is independent of the de¬ 
tails of the crystal structure itself. The eigenvectors cor¬ 
responding to the first few modes will have a structure 
independent of the details of the trap potential, and will 
cause a similar deviation from the power law behavior 
as we saw above. This deviation implies that the spin- 
spin couplings are no longer correlated with the distances 
between the ions, and hence an increase in spatial uni¬ 
formity of the crystal is not guaranteed to have a bear¬ 
ing on the uniformity of the spin-spin interactions. This 
is an important observation, as the increasingly uniform 
nearest-neighbor distances for the triangular wall crystal 
would imply a more uniform spin-spin coupling strength 
only for detuning strength values that are moderately 
large. For intermediate values of S , it is important to 
consider the nature of modes just below the center-of- 
mass mode to describe the spin-spin coupling strength 
between ions corresponding to that strength of detun¬ 
ing. 


IV. CONCLUSIONS AND DISCUSSION 


normalized RMSD = 


max s 


(26) 


as a function of the detuning 6 in Fig. [9] We see an ad¬ 
herence to the power law (characterized by values of the 
normalized RMSD close to 0) for both small and large 
detunings S. More importantly, we see the largest devia¬ 
tion in the normalized RMSD parameter for strengths of 
detuning in the intermediate range of 10 ~ 4 co z to lO 1 ^, 
We can understand this behavior easily if we look at the 
structure of the static part of Eq. that relates the 
strength of the spin-spin coupling, .Ji 3 . to the normal 
mode properties of the crystal. Each term in the sum¬ 
mation can be understood to correspond to an eigen- 
mode’s contribution to the coupling strength. When the 
beat-note frequency /i is very close to to z (6 ~ 0), only 
the center-of-mass mode (corresponding to uniform mo¬ 
tion of the ions) contributes, and the Jij does not de¬ 
pend on distance. This behavior corresponds to a value 
of a = 0. For /i farther away from the center-of-mass 
mode, the lower modes begin to contribute increasingly. 
When only a few modes contribute, we cannot expect 
the power law behavior to hold [TEj . The structure of 
the eigenvectors, as we see in Fig. [Gj is clearly not com¬ 
patible with the power law decay of Jij with distance. 
For many ions sitting at opposite edges, there is a large 
Jij whereas the coupling is virtually zero for other pairs 
separated by much smaller distances. This is the origin 
of the spread of spin-spin couplings in Fig. [7] and the in¬ 
crease in the value of the normalized RMSD in Fig. [9] As 
we move towards larger values of /i (and 5), we see that 
all the modes begin to contribute almost equally, and in 


In this work, we have examined the properties of a 
Penning trap with an additional anharmonic and trian¬ 
gular rotating wall potential which provide a much more 
uniform ionic crystal for use in quantum simulation. By 
performing a detailed analysis of the equilibrium posi¬ 
tions, the phonons, and the effective spin-spin interac¬ 
tions, we find that indeed one can generally obtain more 
uniform spin-spin coupling constants. As one might have 
predicted, the relationship between ionic spacing in the 
lattice and the uniformity of the spin-spin interactions 
is not directly one-to-one. For small a values, it is the 
character of the phonon eigenmodes that lie close to the 
center-of-mass mode that determine the behavior of the 
spin-spin couplings more than the interparticle spacing. 
We hope that the result of this work will be found to 
be useful in planning future experiments with the Pen¬ 
ning trap platform that will employ additional anhar¬ 
monic trap terms and a triangular rotating wall for a 
more uniformly spaced triangular lattice. 
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